p.s <- function(pars =c(0.3,0.2,0.3,0.2,0.1,0.1), time = 1){
	F <- pars[c(1,3)]
	M <- pars[c(2,4)]
	mu <- pars[5:6]

	a <- matrix(c(
		-F[1]-M[1]-mu[1], mu[1], F[1], 0, M[1], 0,
		mu[2], -F[2]-M[2]-mu[2], 0, F[2], 0, M[2],
		rep(0,6*4)), ncol = 6, nrow = 6, byrow = T)
	
	eigs <- eigen(a)
	p <- eigs$vec %*% diag(exp(eigs$val*time)) %*% solve(eigs$vec)
	
	apply(p,1,sum)
	
	a. <- matrix(c(
		-mu[1], mu[1],
		mu[2], -mu[2]), ncol = 2, nrow = 2, byrow = T)
	eigs <- eigen(a.)
	p. <- eigs$vec %*% diag(exp(eigs$val*time)) %*% solve(eigs$vec)
	
	
	return(list(p,p.))
}

p.bad <- function(pars = log(c(0.3,0.2,0.3,0.2,true.p[[2]][1,1]/(1-true.p[[2]][1,1]),
	true.p[[2]][2,2]/(1-true.p[[2]][2,2]))), time = 1){
	e.pars <- exp(pars)
	p.mig <- diag(c(1/(1 + 1/e.pars[5:6]),rep(1,4)))
	p.mig[1,2] <- 1-p.mig[1,1]
	p.mig[2,1] <- 1-p.mig[2,2]
	p.surv <- diag(c(exp(-sum(e.pars[1:2])*time/2),exp(-sum(e.pars[3:4])*time/2), rep(1,4)))
	p.surv[1,c(3,5)] <- e.pars[1:2]*(1-p.surv[1,1])/sum(e.pars[1:2])
	p.surv[2,c(4,6)] <- e.pars[3:4]*(1-p.surv[2,2])/sum(e.pars[3:4])
	out <- p.surv %*% p.mig %*% p.surv %*% p.mig
	return(out)
}
obj.fn <- function(pars, true.pars){
	
	phat <- p.bad(pars, time = 1)
	true <- p.s(true.pars, time = 1)[[1]]
	obj <- sum((log(phat[1:2,]) - log(true[1:2,]))^2)
	return(obj)
}
truth <- c(0.2,0.2,0.2,0.2,0.1,0.1)
true.p <- p.s(truth)
p.bad(log(c(truth[1:4],p.s(truth)[[2]][1,1]/(1-p.s(truth)[[2]][1,1]),
	p.s(truth)[[2]][2,2]/(1-p.s(truth)[[2]][2,2]))))
obj.fn(pars = log(c(0.3,0.2,0.2,0.3,true.p[[2]][1,1]/(1-true.p[[2]][1,1]),
	true.p[[2]][2,2]/(1-true.p[[2]][2,2]))), truth)
x <- optim(par = log(c(truth[1:4],p.s(truth)[[2]][1,1]/(1-p.s(truth)[[2]][1,1]),
	p.s(truth)[[2]][2,2]/(1-p.s(truth)[[2]][2,2]))), obj.fn, true.pars = truth)
truth.range <- expand.grid(seq(0.2,1,0.2),seq(0.2,1,0.2),seq(0.2,1,0.2),seq(0.2,1,0.2),seq(0.2,1,0.2),seq(0.2,1,0.2))
res <- t(apply(truth.range, 1, function(x){
	truth <- x
	p.truth.move <- diag(p.s(truth, time = 1)[[2]])[1:2]
	opt <- optim(par = log(c(truth[1:4],p.truth.move/(1-p.truth.move))), obj.fn, true.pars = truth)
	est <- c(exp(opt$par[1:4]), 1/(1+exp(-opt$par[5:6])))
	true <- c(truth[1:4], p.truth.move)
	out <- c(est, true)
	return(out)}))
apply(abs(res[,7:12]-res[,1:6]),2,max)
par(mfrow = c(2,4))
for(i in 1:6) plot(density(abs(res[,7:12]-res[,1:6])[,i]))
ind <- (abs(res[,7:12]-res[,1:6])[,1]>0.1)
rel.bias <- (res[,7:12]-res[,1:6])/res[,1:6]
library(scatterplot3d)
scatterplot3d(x=truth.range[,1], y = truth.range[,2], z = abs(res[,7:12]-res[,1:6])[,1])
library(rgl)
ind <- which(truth.range[,2]== 0.2 & truth.range[,4] == 0.2) #natmort is 0.2
plot3d(x=truth.range[ind,1], y = truth.range[ind,3], z = rel.bias[ind,1])
plot3d(x=truth.range[ind,1], y = truth.range[ind,5], z = rel.bias[ind,1])
plot3d(x=truth.range[ind,1], y = truth.range[ind,6], z = rel.bias[ind,1])

par(mfrow = c(2,4))
for(i in 1:6) plot(truth.range[ind,i], abs(res[,7:12]-res[,1:6])[ind,1])



bias <- res[,7:12]-res[,1:6])

#for traditional spatial model, suvival and migration are separated temporally:
#P(survive,die) * P(migrate|survive) =

# S1  0 U1  0 M1  0     m11 m12  0  0  0  0    PS11 PS12 PF11 PF12 PM11 PM12
#  0 S2  0 U2  0 M2     m21 m22  0  0  0  0    PS21 PS22 PF21 PF22 PM21 PM22
#  0  0  1  0  0  0       0   0  1  0  0  0       0    0    1    0    0    0
#  0  0  0  1  0  0  x    0   0  0  1  0  0  =    0    0    0    1    0    0 
#  0  0  0  0  1  0       0   0  0  0  1  0       0    0    0    0    1    0
#  0  0  0  0  0  1       0   0  0  0  0  1       0    0    0    0    0    1
#S1 = exp(-(F1 + M1)t)

#In FSCT: P(survive,die,migrate): exp(At) where A =
#  -a1 mu12 F1  0 M1  0
# mu21  -a2  0 F2  0 M2
#    0    0  0  0  0  0
#    0    0  0  0  0  0
#    0    0  0  0  0  0
#    0    0  0  0  0  0
#a1 = M1 + F1 + mu12
# exp(At) = V * exp(U*t) Vinv where U is the vector of eigenvalues and V is the matrix of eigenvectors of A
#set first two rows of exp(At) and P(survive,die) * P(migrate|survive) equal to find bias in F and M estimates
#